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1  Introduction 


Turbulent  flow  is  the  most  frequently  encountered  type  of  viscous  flow,  yet  a  quantitative 
physical  understanding  of  turbulence  remains  an  unsolved  problem.  This  is  true  even  for 
the  simple  case  of  flow  in  a  straight  pipe.  We  do  not  know  how  energy  is  transferred  from 
the  driving  force,  the  pressure  gradient,  to  the  turbulent  fluid  motion.  Hence,  this  simple 
case  provides  a  starting  point  for  investigations  in  turbulent  flow. 

Over  the  years,  many  measurements  have  been  made  in  pipe  flows  and  other  wall- 
bounded  flow  fields,  but  the  conclusions  drawn  from  the  measurements  have  not  been 
final.  While  it  is  possible  to  measure  the  instantaneous  three  components  of  fluid  velocity 
using  laser- Doppler  velocimetry,  this  information  provides  a  very  small  glimpse  of  the 
complete  picture.  Even  with  several  velocity  probes  in  the  flow  field  and  at  the  wall[l], 
the  picture  is  incomplete.  New  experimental  techniques,  such  as  pulsed  laser  velocimetry 
(which  allows  planes  of  instantaneous  velocity  to  be  measured),  are  providing  more  data 
about  the  structure  of  turbulent  wall-bounded  flow  than  any  technique  to  date.  None-the- 
less,  an  accurate  model  of  the  energy  transfer  process  is  still  lacking. 

1.1  Relevance  To  Aerodynamics  and  Ballistics 

Fluid  turbulence  is  important  for  the  computational  fluid  dynamics  solution  of  the  flow  field 
around  a  projectile  in  flight.  Within  this  flow  field  are  wall-bounded  flows,  free  shear  flows, 
and  wake  flows.  If  the  projectile  is  rocket  assisted,  has  a  tracer,  or  contains  a  base-burn 
system,  then  the  flow  field  may  also  contain  jets.  These  four  distinct  types  of  turbulent  flow 
(free  and  wall-bounded  shear  flow,  wakes,  and  jets)  span  the  varieties  of  fluid  turbulence. 
Each  variety  has  turbulence  models,  such  as  mixing  length  or  production-dissipation,  that 
are  based  on  empirical  data  and  require  tuning  for  a  specific  application.  Direct  simulation 
has  no  such  models,  but  has  severe  limits  in  the  complexity  of  the  geometry,  the  size  of 
the  Reynolds  numbers,  and  the  time  required  to  achieve  fully  developed  turbulent  flow. 

In  spite  of  the  limits,  direct  simulation  provides  information  that  is  unavailable  from 
physical  experiments,  namely  the  full  three-dimensional  time-varying  flow  field.  In  this 
way,  direct  simulation  is  a  “computer  experiment"  that  can  be  used  to  tune  turbulence 
models.  The  long  range  goal,  however,  is  to  develop  more  complete  models  for  fluid  tur¬ 
bulence  needed  by  engineers  to  compute  the  mean  flow  field  surrounding  a  projectile.  The 
steps  to  this  goal  include  using  the  direct  simulation  to  obtain  a  flow  field  for  comparison 
to  large  eddy  simulation  of  the  small  scale  motion.  This  may  lead  to  better  understanding 
and,  ultimately,  to  better  models  of  fluid  turbulence. 

1.2  The  Universal  Law  of  the  Wall 

The  early  measurements  of  velocity  in  wall-bounded  flows  allowed  the  construction  of  the 
“Universal  Law  of  the  Wall”  that  is  depicted  in  Figure  1.  This  figure  shows  that  the 
streamwise  fluid  velocity  component  is  zero-valued  at  the  wall.  W'ithin  the  “viscous  wall 
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Figure  1:  Regions  in  Wall  Bounded  Flows 

layer”  (so  called  because  of  the  strong  effects  of  viscosity),  the  velocity  increases  to  a  value 
near  that  of  the  time-averaged  maximum  fluid  velocity.  The  region  beyond  the  viscous 
wall  layer  is  referred  to  as  the  outer  flow.  As  defined  in  this  report,  the  viscous  wall  layer 
contains  the  viscous  (also  called  laminar  or  linear)  sublayer  and  the  buffer  layer,  which 
are  commonly  mentioned  in  the  literature[2].  The  outer  flow  is  sometimes  called  the  core 
region,  and  it  is  a  combination  of  the  inertial  (or  logarithmic)  sublayer,  and,  if  present, 
the  wake  region  [2]. 

At  first,  researchers  thought  that  the  instability  that  made  a  flow  become  turbulent  was 
to  be  found  mainly  in  the  outer  flow.  They  suggested  that  the  velocity  fluctuations  found 
in  the  outer  flow  were  associated  with  large  eddy  structures.  These,  in  turn,  fed  energy 
to  smaller  and  smaller  eddy  structures  until,  at  the  wadi,  the  smallest  eddy  structures 
transferred  the  energy  to  the  wall  as  heat.  In  this  sense,  the  outer  flow  controlled  the 
motion  of  the  fluid  in  the  viscous  wall  layer.  The  fluid  near  the  wall  responded  to  the 
velocity  fluctuations  of  the  outer  flow  without  influencing  the  outer  flow.  This  model  is 
not  “wrong;”  the  large  eddies  do  feed  the  smaller  eddies.  However,  later  research  showed 
that  the  motion  of  the  fluid  in  the  viscous  wall  layer  might  have  a  strong  influence  on  the 
fluid  in  the  outer  flow. 
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1.3  Turbulent  Kinetic  Energy 

At  present,  the  model  of  fluid  turbulence  in  a  wall  bounded  flow  surrounds  the  two  main 
contributors  to  the  turbulent  kinetic  energy  budget  (see  Section  6) — the  production  and 
dissipation  of  turbulent  kinetic  energy.  Researchers  have  found  that  the  production  and 
dissipation  do  not  add  to  zero  in  the  viscous  wadi  layer;  there  is  extra  turbulent  kinetic 
energy  available.  The  result  is  a  supply  of  turbulent  energy  to  the  outer  flow,  where  it  is 
dissipated.  Thus,  the  velocity  fluctuations  in  the  outer  flow  can  result  from  an  influx  of 
turbulent  kinetic  energy  from  the  viscous  wall  layer.  This  is  a  fax  different  notion  than  the 
earlier  model  which  suggested  that  the  fluid  in  the  viscous  wall  layer  merely  responded  to 
the  fluctuations  of  velocity  in  the  outer  flow.  Recent  measurements  show  that  the  viscous 
wall  layer  contributes  to  those  fluctuations;  in  fact,  the  motion  of  the  fluid  in  the  viscous 
wall  layer  may  control  the  motion  of  the  outer  flow.  The  eddies  in  the  viscous  wall  layer 
still  seem  to  be  influenced  by  the  outer  flow.  This  suggests  that  there  may  be  a  feedback 
involved  in  the  process.  This  model,  which  is  far  from  the  final  model,  is  the  starting  point 
for  the  present  work. 

Currently,  the  interest  is  in  developing  techniques  for  examining  the  eddies  in  the 
viscous  wall  layer.  Instead  of  a  pipe  flow,  however,  the  present  study  is  performed  on  the 
flow  field  between  two  parallel,  flat  plates.  Thus,  the  flow  field  geometry  is  more  closely 
related  to  channel  flow,  although  the  features  of  the  flow  field  are  quite  similar  to  pipe 
flow.  This  paper  documents  the  initial  verification  that  the  simulation  reproduces  nature. 

1.4  The  Present  Study 

As  expected  from  the  periodicity[3]  in  turbulent  flow,  Fourier  methods  arise  naturally  in 
the  analysis  of  turbulence  and  they  axe  exploited  in  the  present  study.  A  direct  numerical 
simulation  of  turbulent  channel  flow,  using  spectral  techniques,  has  been  developed[4],  and 
the  results  are  compared  to  experimental  measurements.  Since  many  models  exist  that 
predict  mean  quantities  in  turbulent  channel  flow,  the  most  accurate  of  those  models  is 
also  discussed  as  a  comparison  to  the  direct  simulation. 

Research  on  turbulent  fluid  flow  has  been  invigorated  in  the  past  ten  years  by  the  devel¬ 
opment  of  supercomputer  simulations  of  flow  fields,  although  direct  numerical  simulation 
of  turbulent  flow  in  a  simple  geometry  is  still  a  formidable  task.  The  difficulty  resides  in 
the  wide  range  of  spatial  scales;  hence,  the  computational  domain  must  be  sized  to  capture 
the  largest  scale  of  motion  while  at  the  same  time  have  small  enough  grid  spacing  to  resolve 
the  smallest  relevant  scale.  Three  dimensional  effects  must  also  be  considered  because  of 
the  nonlinear  terms  in  the  Navier-Stokes  equations,  v*Vv.  These  represent  the  conversion 
of  linear  momentum  into  angular  momentum,  or  the  formation  of  vorticity.  To  account 
for  vorticity,  calculations  in  three  dimensions  are  required,  even  for  channels  where  the 
time-averaged  flow  is  two-dimensional. 
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2  Background 


2.1  Reynolds  Decomposition  of  the  Navier-Stokes  Equations 

The  equations  governing  the  motion  of  an  incompressible  fluid  are  the  Navier-Stokes  equa¬ 
tions: 

^  =  — vVv--Vp  +  i/V2v  (1) 

at  p 

and  the  continuity  equation: 

V-v  =  0.  (2) 

In  these  equations,  v  is  the  instantaneous  three  dimensional  velocity  vector,  t  is  the  time, 
p  is  the  density  of  the  fluid,  p  is  the  hydrodynamic  pressure,  and  v  is  the  kinematic  viscos¬ 
ity.  Here  the  coordinate  system  is  Cartesian,  and  (xj, x2,x3)  are  the  streamwise.  normal, 
and  spanwise  directions  in  the  channel  geometry  with  corresponding  velocity  components 
(vi,v2,v3). 

If  the  instantaneous  velocity  vector  is  replaced  in  Equations  1  and  2  with  a  velocity 
vector  that  is  decomposed  into  a  mean  flow,  V,  and  a  fluctuating  velocity,  v,  such  that: 


v  =  V  +  v 


(3) 


then  the  resulting  “Reynolds-decomposed”  Navier-Stokes  equations  contain  an  additional 
tensor: 

t,j  =  - pvii’j •  (4) 

Here,  rt]  is  the  Reynolds  stress  tensor,  v,  and  f;  are  the  components  of  the  instantaneous 
velocities,  and  the  overbar  represents  the  time  average. 

This  decomposition  into  mean  flow  and  turbulent  velocity  fluctuations  isolates  the 
effects  of  the  fluctuations  from  the  mean  flow.  Now,  however,  there  are  nine  components 
of  Tij  (six  are  independent)  in  addition  to  the  p  and  the  three  components  of  V.  Herein 
lies  the  closure  problem  of  turbulence.  While  this  produces  additional  stress  terms,  it  does 
allow  solution  of  the  time-averaged  Navier-Stokes  equations. 

One  common  method  of  closing  the  equations  is  to  use  the  similarity  of  the  form  of  the 
Reynolds  stress  to  the  viscous  stress  in  the  Navier-Stokes  equations.  An  eddy  viscosity,  vt, 
is  defined  by: 

dU{x 2)  ,  s 

r12  =  pvt—r -  (5) 

OX  2 

where  U(x2 )  is  the  streamwise  component  of  the  mean  velocity. 


2.2  Channel  Flow 

The  best  empirical  fit  of  the  velocity  profile  in  a  flat-walled  channel  can  be  represented  by 
the  function  for  the  eddy  viscosity  developed  by  Reynolds  and  Tiederman[5): 

\  +  ~  =  +  F2G2  (6) 

m  4 
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where 


and 


In  this  equation  h  is  the  half-height  of  the  channel,  r  is  the  local  shear  stress,  k  is  the  von 
Karman  parameter,  A o  is  the  van  Driest  parameter,  and  tw  is  the  shear  stress  at  the  wall 
defined  as: 

_  dU(x  2) 


Tw  =  pU- 


dx<2 


(7) 


X2=0 


When  this  value  for  the  eddy  viscosity  is  substituted  into  the  mixing  length  model  for 
the  Reynolds  stress,  and  that  in  turn  is  used  in  the  Reynolds-decomposed  Navier-Stokes 
equations,  an  excellent  approximation  to  the  time- averaged  flat- walled  channel  velocity 
profile  results. 

Unfortunately,  this  result  tells  us  nothing  about  the  dynamics  of  the  flow  field;  specif¬ 
ically,  it  does  not  describe  the  instantaneous  interaction  between  the  wall  and  the  outer 
flow.  Since  it  is  based  on  empirical  data,  an  experiment  is  generally  required  to  tune 
the  parameters[6,7,8].  Direct  numerical  simulation  has  no  such  requirement  and  has  the 
advantage  of  generating  instantaneous  velocities  in  three  dimensions.  This  produces  a 
“computer  experiment,”  that  can  focus  other  experimental  techniques,  such  as  pulsed 
laser  velocimetry[9]  and  multi-probe  hot  wire  anemometryf  1,1 0,11].  At  present,  direct 
numerical  simulation  is  limited  to  simple  geometries  and  low  Reynolds  numbers. 


2.3  Viscous  Units 


It  is  convenient  to  mention  viscous  units  here,  which  result  from  making  velocity  and 
hydrodynamic  pressure  dimensionless.  This  is  accomplished  with  the  wall  parameters 
kinematic  viscosity,  u,  and  skin-friction  velocity,  v*.  Viscous  units  are  denoted  by  giving 
the  quantity  a  plus  sign  superscript. 

Thus,  hydrodynamic  pressure,  velocity,  and  length  will  be  made  dimensionless  in  the 
following  manner: 


(8) 


This  allows  comparison  with  experimental  results  that  do  not  have  the  same  half-height 
as  the  channel  in  the  computer  experiment. 


3  A  Model  of  Wall  Turbulence 

Although  the  velocity  fluctuations  in  a  turbulent  flow  field  are  random,  coherent  structures 
exist  in  the  viscous  wall  layer.  This  was  first  observed  by  Beatty,  Ferrel  and  Richardson[12] 
who  showed  “streaky  structures”  that  were  oriented  with  the  mean  flow  and  repeated 
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Figure  2:  Conceptualization  of  Coupled  Near- Wall  Structures 

periodically  in  the  direction  perpendicular  to  the  mean  flow  (the  spanwise  direction  in 
channel  flow).  These  streaky  structures,  marked  with  dye  or  hydrogen  bubbles,  were 
found  to  be  regions  of  low  velocity  fluid;  hence,  the  fluid  had  low  momentum.  The  region 
between  the  streaky  structures  contained  high  momentum  fluid. 

Streaky  structures,  marked  with  dye,  traveled  downstream  and  then  were  lifted  away 
from  the  viscous  wall  layer  toward  the  outer  flow  where  they  dissipated.  This  process  has 
been  referred  to  as  bursting,  and  it  occurs  periodically. 

Quantifying  bursting  and  the  spacing  of  streaky  structures  has  been  accomplished  by 
many  researchers  (see  Blackwelder  and  Haritonidis[13]).  These  quantities  are  independent 
of  Reynolds  number  and  are  found,  when  made  dimensionless  with  wall  parameters,  to 
be  A+  ss  100  for  the  spanwise  spacing  of  the  streaky  structures.  The  length  of  a  streaky 
structure  is  roughly  10  times  its  width[14].  Further  experimental  work  showed  that  the 
streaks  were  associated  with  eddies. 

In  an  effort  to  explain  this  behavior,  Hanratty[14]  developed  the  picture  first  presented 
by  Townsend[15].  Hanratty  suggests  that  these  eddies  have  inflows  and  outflows  that  are 
of  the  same  strength  and  that,  in  a  time-averaged  sense,  the  in-  and  outflows  are  coupled, 
as  shown  in  Figure  2.  The  eddies  create  streamwise  velocity  fluctuations  by  bursting  and 
bringing  low  momentum  fluid  from  the  wall  into  the  outer  flow.  The  opposite  is  also  true; 
i.e.,  they  carry  high  momentum  fluid  toward  the  wall.  High  momentum  fluid  is  transported 
by  negative  (toward  the  wall)  normal  velocities;  low  momentum  fluid  is  transported  by 
positive  normal  velocities.  Hence,  with  this  transfer  process,  when  the  normal  velocity  is 
positive,  the  fluctuating  streamwise  velocity  is  negative.  This  means  that  both  in-  and 
outflows  are  associated  with  large  values  of  negative  Reynolds  shear  stress,  —viV2- 
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4  Governing  Equations  and  Method  of  Solution 

The  Navier-Stokes  equations  and  the  continuity  equation  (Equations  1  and  2)  can  be 
rewritten  using  the  following  definitions: 

<*>  =  V  x  v  ;  ITr  =  —  +  -  +  ;V2;  v 2  =  v-v  (9) 

P  P  2 

where  u  is  the  vorticity.  Term  IIj  is  the  total  pressure  consisting  of  the  mean,  P,  and 
fluctuating,  p,  hydrodynamic  pressure  and  the  dynamic  pressure.  The  advection  term  in 
the  Navier-Stokes  equations  becomes: 

v*Vv  =  —v  x  u>  +  ^Vv2  (10) 

and  the  Navier-Stokes  equations  can  be  written  with  vorticity  explicit: 

C^V 

—  =  v  x  u>  —  vrij  +  vV2\.  (ii) 

at 

This  form  of  the  equations  makes  the  numerical  solution  somewhat  easier. 

4.1  The  First  Fractional  Time  Step 

A  three-level  fractional  step  method,  developed  by  Orszag[l6],  is  used  to  solve  the  equa¬ 
tions.  This  amounts  to  setting  the  left-hand  side  of  Equation  11  equal  to  each  one  of  the 
terms  on  the  right-hand  side  separately  (hence,  three  steps). 

By  subtracting  the  mean  streamwise  hydrodynamic  pressure  gradient  (the  mean  span- 
wise  and  normal  gradients  are  zero)  from  the  second  term  on  the  right-hand  side  of  Equa¬ 
tion  11  and  modeling  the  streamwise  advective  term  implicitly  to  reduce  convective  insta¬ 
bility,  the  advection  part  of  Equation  11  becomes: 

d\  dv  dv 

|  +  t/(X!)_  =  W  +  t:(l2)_,  02) 

where  U(z :?)  is  the  streamwise  component  of  the  mean  velocity.  Function  H  is  equal  to 
v  x  u>  —  ,  where  ~  is  the  mean  streamwise  hydrodynamic  pressure  gradient.  Solving 
Equation  12  by  application  of  a  Crank-Nicholson  scheme  to  the  left  side  and  an  Adams- 
Bashforth  scheme  to  the  right  side  gives  the  following  result: 


.  W) {?fl 

At  2  2  2  dx\  dx\  dx\ 


the  superscript  refers  to  the  time  step  and  the  hat  (e.g.,  v)  represents  quantities  resulting 
after  the  first  fractional  time  step,  also  called  the  “inviscid  step."  Solving  for  v"+1  gives 
the  velocity  after  the  first  fractional  time  step. 


4.2  The  Second  Fractional  Time  Step 

The  no-penetration  boundary  condition,  v?  =  0  at  the  wails,  is  applied  at  the  second 
fractional  step.  An  implicit  Euler  scheme  is  used  which  results  in: 

£n+i  _  yn+1  =  -At  Vrin+1.  (14) 


Hence  the  second  step  imposes  incompressibility  and  the  effects  of  pressure,  and  the  double 
hat  refers  to  the  velocity  after  this  step.  The  total  pressure  (Hx)  was  modified  for  the  first 
step  by  subtracting  the  mean  hydrodynamic  pressure.  Now  term  II  contains  the  fluctuating 
hydrodynamic  pressure,  p,  and  the  dynamic  pressure  so  the  subscript  T  will  no  longer  be 
used: 

n  =  £  +  ^vJ.  (is) 

P  2 

The  original  method  of  Orszag’s  did  not  properly  enforce  the  boundary  condition  for 
the  pressure  field,  and  continuity  was  not  satisfied  at  the  walls.  Using  a  Green’s  function 
method  (see  Lyons  [4])  in  this  second  step  corrects  the  calculated  pressure  field  and  satisfies 
continuity.  The  pressure  term  is  determined  by  requiring  the  divergence  of  the  velocity 
field  to  be  zero  after  the  second  fractional  step: 


Thus,  from  14: 

with  the  boundary  condition: 


v2rin+1  = 


dnn+1 


dxi 


=  v  ■ 


rj  =0 


=  0. 

(16) 

V*vn+1 

(17) 

At 

d2vn+1.j1 

(IB) 

dx\ 

This  is  a  linear  inhomogeneous  equation  with  inhomogeneous  boundary  conditions.  The 
equation  can  be  split  into  two  parts  (an  inhomogeneous  equation  with  homogeneous  bound¬ 
ary  conditions  and  a  homogeneous  equation  with  inhomogeneous  boundary  conditions). 
Additionally,  the  pressure  (II)  is  expressed  as  the  sum  of  an  inviscid  and  a  viscous  pres¬ 
sure.  The  viscous  part  of  the  pressure  satisfies  the  homogeneous  equation,  and  its  boundary 
conditions  are  satisfied  if  and  only  if  continuity  is  satisfied  at  the  walls.  In  other  words, 
the  inhomogeneous  boundary  condition  is  satisfied  if  and  only  if  Equation  16  is  true  at 
the  walls.  The  inviscid  part  satisfies  the  inhomogeneous  equation,  although  this  does  not 
insure  continuity  is  satisfied. 


4.3  The  Third  Fractional  Time  Step 

The  final  step  takes  viscosity  into  account  and  imposes  the  no-slip  boundary  conditions, 
vn+1  =  0,  at  the  walls.  It  is  an  implicit  backward  Euler  scheme: 

v»+>  _*«-H  =uAt-  V2vn+1.  (19) 
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Non-wall  boundary  conditions  are  periodic;  the  outflow  at  one  “end”  is  introduced  as  the 
inflow  at  the  other  “end.”  This  assumption  makes  the  problem  suitable  for  a  Fourier  or 
“spectral”  approach.  (It  is  reasonable  if  the  dimensions  in  the  spanwise  and  streamwise 
computational  domain  are  more  than  twice  the  largest  relevant  scale  of  motion  in  these 
directions.  This  information  is  obtained  by  experiment,  and  is  discussed  below.) 


4.4  The  Temperature  Field 

Additionally,  a  temperature  gradient  is  imposed  on  the  flow  field.  The  walls  are  held  at 
constant  temperatures  that  are  different  from  each  other.  The  temperature  field  is  then 
solved  in  the  first  and  third  fractional  steps  of  the  computation.  The  temperature  equation 

is  ar 

_  =  _v.VT  +  —  V2T  (20) 

where  T  is  the  temperature  and  Pr  is  the  Prandlt  number,  which  is  a  ratio  of  the  kinematic 
viscosity  (momentum  diffusivitv)  to  the  thermal  diffusivity.  The  thermal  diffusivity  is  the 
ratio  of  the  thermal  conductivity  to  the  density  and  heat  capacity. 

In  order  to  get  homogeneous,  Dirichlet  boundary  conditions  the  temperature  is  trans¬ 
formed  by 

T  =  6-Ax2  (21) 

where  6  is  the  transformed  temperature  and  A  is  a  constant  that  controls  the  temperature 
difference  between  the  two  walls.  Following  the  same  procedure  as  for  the  convective  (first) 
step  of  the  Navier-Stokes  solution  scheme,  and  denoting  —  v  •  V9  —  Av2  as  I,  the  result  is 


Qn+\  _  Qn 

At 


-  1 


d0n+ 1 

dx\ 


o?t 

~  di\ 


+ 


den~x 


dx\ 


(22) 


The  solution  scheme  is  the  same  as  for  Equation  13. 

Since  there  is  no  pressure  term  in  the  temperature  equation,  the  remaining  terms  are 
solved  in  the  third  fractional  step  as  in  Equation  19: 

0n+1  _  fln+i  _  .  V20n+1,  (23) 

Pr 

and  the  boundary  conditions  axe  homogeneous. 


4.5  The  Spectral  Technique 

Using  the  spectral  method  of  Orszag[16],  the  velocity  field,  v  =  (iq ,  v2,  t’3).  can  be  trans¬ 
formed  using  the  following  expansion: 

N 

v(x,t)  =  H  £  rim„(f )  exp 

l'l<£  M<f  n=0 
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x2 

k‘ 


Figure  3:  Channel  Geometry 

In  this  equation,  L  and  M  are  the  number  of  Fourier  modes  in  the  jj  and  13  directions,  Tn 
is  a  Chebyshev  polynomial  in  the  xi  direction  and  N  is  the  degree  of  the  polynomial.  Terms 
Ai  and  A3  are  the  lengths  of  the  computational  domain  in  the  streamwise  and  spanwise 
directions,  and  h  is  the  channel  half-height  as  shown  in  Figure  3.  The  Fourier  coefficients 
of  the  expansion  are  represented  by  t>(mn(0<  and  these  are  referred  to  collectively  as  the 
velocity  field  in  spectral  space.  The  hydrodynamic  pressure  and  the  temperature  fields  are 
expanded  in  a  similar  manner. 

This  transformation  reduces  the  problem  from  solving  partial  differential  equations  for 
the  velocity,  pressure,  and  temperature  fields  to  solving  equations  for  the  spectral  coeffi¬ 
cients  of  the  physical  fields.  Spatial  derivatives  of  the  velocity  field  are  easily  calculated 
in  spectral  space;  the  first  derivative  of  the  velocity  field  in  spectral  space  with  respect  to 
the  streamwise  direction  is  equal  to  the  the  original  coefficients  multiplied  by  the  factor 
2nil/Xx.  Derivatives  with  respect  to  the  normal  direction  are  slightly  more  complicated, 
but  do  not  require  finite  differencing  and  are  exact. 

The  coefficients  of  the  Fourier  and  Chebyshev  series  expansions  are  determined  by 
forcing  the  series  to  be  an  exact  representation  of  the  corresponding  velocity,  temperature, 
or  pressure  field  at  the  grid  points.  A  Chebyshev  series  can  be  rewritten  as  a  Fourier 
series;  this  is  used  to  re-arrange  Equation  24  to  be  the  sum  of  three  Fourier  series.  This 
allows  calculation  of  the  physical  space  fields  by  fast  Fourier  transform — and  thus  allows 
calculation  of  the  spectral  space  coefficients  by  inverse  fast  Fourier  transform. 

4.6  The  Solution  Scheme 

In  the  solution  scheme,  the  spectral  representations  of  the  velocity,  temperature,  and 
pressure  fields  are  substituted  into  the  various  finite  difference  equations  in  the  three 
fractional  steps. 
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In  the  first  fractional  step,  the  terms  in  function  H  are  computed  in  physical  space. 
This  is  done  by  transforming  the  velocity  and  vorticity  fields  to  physical  space  via  fast 
Fourier  transform,  then  transforming  the  resulting  “if”  field  back  to  spectral  space.  Also, 
this  procedure  is  used  for  the  first  fractional  time  step  in  the  temperature  field.  These  terms 
and  the  spectral  expansions  for  the  velocity  and  temperature  fields  axe  then  substituted 
into  Equations  13  and  22. 

The  second  and  third  fractional  steps  do  not  require  terms  to  be  calculated  in  physical 
space.  In  these  steps,  the  spectral  expansion  is  substituted  into  the  appropriate  finite 
difference  equation  and,  with  some  manipulation,  the  values  of  the  various  fields  are  cal¬ 
culated  at  the  next  time  step.  Since  the  Navier-Stokes  equations  describe  the  motion  of 
the  fluid,  the  dynamics  of  the  flow  field  are  reproduced  numerically. 

As  described,  this  spectral  method  produces  a  finite  grid  of  points  with  sine  waves, 
representing  the  velocity,  superimposed  upon  it.  These  sine  waves  are  multiplied  together 
in  the  vorticity-velocity  cross  product  of  the  Navier-Stokes  equations,  resulting  in  sums 
and  differences  of  wavenumbers  (inverse  wavelength).  Longer  wavelengths  can  be  resolved 
on  the  grid,  but  the  shorter  wavelengths  can  not  be  resolved  and  are  aliased  onto  longer 
wavelengths.  (Note  that  the  summations  in  Equation  24  are  centered  on  zero;  hence, 
infinite  wavelength  can  be  resolved  but  wavenumbers  greater  than  half  the  number  of 
Fourier  modes  can  not  be  resolved.)  In  nature,  the  energy  of  the  short  wavelengths  would 
produce  heat,  not  kinetic  energy.  To  rid  the  simulation  of  the  short  wavelengths,  the  top 
third  of  the  wavenumbers  in  the  streamwise  and  spanwise  directions  are  set  to  zero  before 
forming  the  products.  This  is  referred  to  as  “dealiasing.”  It  would  appear  that  half  of  the 
wavenumbers  should  be  set  to  zero  in  order  to  dealias,  but  Orszag[17]  has  shown  that  one 
third  is  sufficient. 

4.7  Other  Research 

Two  other  research  groups  have  worked  on  direct  numerical  simulation  and  have  developed 
slightly  different  solution  schemes.  Moin  and  Kim[18]  typically  have  used  finite  difference 
schemes  in  place  of  the  Chebyshev  polynomial.  Similarly,  Kleiser  and  Schumann[19]  use 
finite  differences,  although  in  the  solution  they  employ  sequences  of  one-dimensional  scalar 
Helmholtz  equations  instead  of  solving  block-tridiagonal  matrices. 

To  date,  the  results  of  Moin  and  Kim[18]  and  Orszag  and  Kells[17]  appear  to  reflect 
turbulent  channel  measurements  the  best.  Orszag  and  Kells  were  more  interested  in  the 
transition  from  laminar  to  turbulent  flow;  consequently,  their  algorithm  is  not  suitable  for 
higher  Reynolds  numbers.  The  work  of  Moin  and  Kim  has  approximately  the  resolution 
of  this  current  study. 

5  Computational  Parameters 

In  this  numerical  study,  the  dimensionless  half-height  of  the  channel  was  150  viscous  units 
for  a  Reynolds  number  of  2260  based  on  this  characteristic  length  scale  (approximately  a 
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pipe  Reynolds  number  of  8,000).  The  Reynolds  number  is  ^  where  Vb  is  the  bulk  velocity, 
or  the  integral  of  the  average  velocity  over  the  height  of  the  channel.  To  resolve  the  large- 
scale  eddies,  the  spanwise  and  streamwise  dimensions  were  set  at  950  and  1900  viscous 
units.  This  produced  a  bulk  velocity  of  13.8  in  viscous  units.  The  number  of  spectral 
modes  for  the  results  shown  in  this  paper  were  (L,  JV,  M)  =  (128, 128, 128). 

The  Reynolds  number  determines  the  largest  and  smallest  scales  of  motion  in  the  flow 
field,  and  it  is  fixed  by  the  half- height,  h+.  This  can  be  shown  by  a  force  balance  over 
the  computational  domain  (see  Figure  3).  Since,  in  the  spanwise  direction,  the  average 
velocity  is  zero  and  there  is  no  hydrodynamic  pressure  gradient,  the  force  balance  is: 

2h\3(P\l  -  P\3)  +  pfcA3(^|i  -  ^|2)  -  2Aj A3  \tv\  =  0.  (25) 

This  can  be  reduced  because  the  average  velocity  is  the  same  at  cross  sections  1  and  2, 
and  the  hydrodynamic  pressure  gradient  is  constant: 


dP+  _  J_ 
dx i+  h+ ' 


(26) 


Thus  picking  a  half- height  determines  the  mean  pressure  gradient,  which  is  the  driving 
force.  The  pressure  gradient  fixes  the  Reynolds  number  since  average  velocity,  kinematic 
viscosity,  and  the  length  scale  are  fixed. 

At  this  point,  the  other  lengths  in  the  computational  domain  may  be  determined.  A 
periodic  boundary  condition  can  be  used  if  the  computational  domain  is  at  least  twice  the 
size  of  the  largest  periodic  structure.  The  appropriate  functions  for  determining  this  size 
are  the  two-point  velocity  correlation  functions,  one  in  the  streamwise  direction: 


f?,,(Axi,0,0)  =  vi(xi,x2,x3)vi(x1  +  Axi,x2,x3) 


(27) 


and  a  similar  function  in  the  spanwise  direction.  From  experimental  data[20]  the  two-point 
velocity  correlation  functions  vanish  at  3.2 h+  in  the  streamwise  direction  and  1.6h+  in  the 
spanwise  direction,  while  others[l8]  suggest  values  of  6h+  and  3 h+.  Before  vanishing,  these 
two-point  velocity  correlation  functions  become  negative,  and  at  this  spanwise  distance 
distinct  streamwise  streaks  form  near  the  wall  in  hydrogen  bubble  experiments [3]. 

The  small  scales  of  motion  must  also  be  considered  in  order  to  determine  the  minimum 
number  of  grid  points  or  Fourier  modes.  Relevant  scales  would  begin  with  the  Taylor 
micro-scale,  which  fixes  the  rate  of  dissipation  of  turbulent  kinetic  energy,  and  increase  in 
size.  The  smallest  scales,  the  Kolmogoroff  scales,  are  purely  dissipative  and  thus  are  not 
important  in  the  direct  simulation,  although  they  are  crudely  represented  by  dealiasing. 
Experimental  measurements[ll]  show  the  Taylor  micro-scale  to  be  about  12  viscous  units 
in  the  spanwise  direction  and  ten  times  that  in  the  streamwise  direction  at  the  wall,  where 
the  scales  are  smallest.  Table  1  shows  that  the  Taylor  micro-scale  may  not  be  resolved  at 
the  wall  in  the  spanwise  direction  for  this  calculation.  Results  suggest  that  this  is  a  minor 
shortcoming. 

Using  these  values,  the  size  of  the  computational  domain  and  the  grid  spacing  can  be 
determined.  Table  1  shows  the  resolution  available  in  the  current  study.  The  modes  were 
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Direction 

Length 

Modes 

Resolution 

Resolvable  Scales 

Spanwise 

950 

85 

11.18 

22.4  to  475 

Stream  wise 

1900 

85 

22.35 

44.7  to  950 

Normal 

300 

128 

3.68 

.01  to  300 

|  Length,  resolution  and  resolvable  scales  in  viscous  units. 

Table  1:  Computational  Domain  Resolution 


adjusted  for  dealiasing.  Also,  in  the  normal  direction  the  spacing  is  stretched  toward  the 
channel  center  so  that  near  the  wall,  where  the  smallest  scales  are  present,  the  resolution 
is  highest. 


6  Method  of  Analysis 

6.1  General  Averages  of  Turbulence  Quantities 

The  solution  procedure  outlined  above  generates  a  computer  experiment  that  is  similar 
to  physical  experimented  facilities,  such  as  water  tunnels  using  particle  image  velocimetry 
and  laser-Doppler  velocimetry.  One  advantage  of  the  computer  experiment  is  the  relative 
simplicity  of  obtaining  instantaneous,  three  dimensional  velocity  measurements  at  any 
number  of  locations,  simultaneously,  in  the  flow  field. 

In  the  present  work  the  reduction  of  data  that  would  accompany  a  physical  experiment 
has  been  performed,  including  averages  of  velocities  and  moments  of  velocity.  For  a  physical 
experiment  the  averages  would  be  taken  either  as  an  ensemble  or  in  time.  Thus,  at  one 
point  in  the  flow  field  the  mean  velocity  would  be  given  by: 

1  rl  o+1' 

V  =  lim  —  /  \  dt  (28) 

‘c-»oo  tt  Jt0 

where  tt  is  some  elapsed  time.  A  spatial  average  is  used  in  this  work  such  that  a  plane  of 
grid  points  in  physical  space,  parallel  to  the  solid  walls,  forms  the  basis  of  the  average.  To 
get  higher-order  statistics,  several  planes  are  averaged  together. 

Averages  of  the  fluctuating  velocities  and  their  moments  will  be  presented.  The  normal 
and  spanwise  mean  velocities  are  zero.  The  second-order  moments  are  the  components  of 
the  Reynolds  stress,  but  are  also  called  turbulence  intensities.  The  third-  and  fourth-order 
moments  are  referred  to  as  the  skewness  (viiijVk)  and  flatness  (v*v}vk).  Skewness  and 
flatness,  involving  velocities  in  two  directions,  are  not  well  understood,  and  therefore  only 
v?  and  v *  are  presented. 
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6.2  The  Turbulent  Energy  Budget 


Turbulent  energy  production  and  dissipation  are  also  of  interest.  An  equation,  commonly 
called  the  turbulent  kinetic  energy  budget,  that  governs  the  mean  kinetic  energy,  of 

the  turbulent  velocity  fluctuations  can  be  obtained  from  the  Navier-Stokes  equations.  Mul¬ 
tiplying  Equation  1  by  the  fluctuating  velocity,  taking  the  time  average,  and  subtracting 
the  mean  kinetic  energy  budget  gives: 


l_ 


-ViVjSij  -  —  ^ VjP  -  2 uvtSij  )  -  2 wSijSij. 


(29) 


The  mean  strain  rate,  Sij ,  and  the  strain- rate  fluctuations,  itJ ,  are  the  symmetric  parts  of 
the  mean  and  fluctuating  deformation  rates.  The  total  deformation  rate  is 


di\  _  _  1  /  dv,  duA  1  ( dv,  dt>A 

dij  S':  r':  ~  2  dx, )  +  2  dxi ) 


(30) 


and  the  mean  and  fluctuating  rates  are  similarly  defined. 

Equation  29  shows  that  the  rate  of  change  of  the  turbulent  kinetic  energy  is  due  to  two 
kinds  of  deformation  work  and  three  forms  of  transport.  If  the  energy  flux  in  a  control 
volume  is  zero,  two  of  the  transport  terms  (the  second  and  third  terms  on  the  right-hand 
side)  redistribute  the  energy  from  one  location  to  another.  The  pressure  gradient -velocity 
interaction  (the  first  term  on  the  right-hand  side)  arises  from  pressure  fluctuations.  In 
an  incompressible  fluid  the  contribution  of  the  mean  pressure  to  the  deformation  work  is 
zero,  but  the  fluctuating  pressure  contributes  to  the  turbulent  kinetic  energy.  Thus,  the 
first  two  terms  represent  the  interaction  of  the  velocity  with  the  sum  of  the  fluctuating 
hydrodynamic  pressure  and  the  dynamic  pressure. 

Deformation-work  terms,  the  last  two  terms  on  the  right-hand  side,  represent  turbu¬ 
lence  production  and  viscous  dissipation.  Turbulence  production  describes  the  exchange 
of  kinetic  energy  between  the  mean  flow  and  the  fluctuations.  Viscous  dissipation  is  the 
rate  at  which  viscous  stresses  perform  deformation  work  against  the  fluctuating  strain  rate. 
These  terms  are  of  interest  in  an  instantaneous  sense,  too,  as  described  in  the  introduction. 

For  plane  channel  flow,  Equation  29  may  be  written  as: 


0  = 


(31) 


The  five  terms  on  the  right  side  of  this  equation  are  the  production,  turbulent  transport, 
fluctuating  pressure  gradient-velocity  interaction  transport,  rate  of  viscous  diffusion,  and 
the  isotropic  dissipation  of  turbulent  kinetic  energy,  respectively.  They  are  hereafter  de¬ 
noted  as  V,T,Q,V,  and  V ,  respectively.  In  this  more  specific  case,  terms  T  and  Q  are 
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associated  with  only  the  normal  velocity  fluctuations,  and  term  V  is  affected  only  by  the 
change  in  the  normal  direction.  Term  V  is  reduced  to  the  Reynolds  shear  stress  interaction 
with  the  mean  flow  velocity  gradient,  and  this  is  the  major  source  for  turbulent  kinetic 
energy.  This  shows  the  importance  of  the  Reynolds  shear  stress.  Since  term  V  is  isotropic 
(in  its  expansion,  there  are  nine  terms),  it  remains  unchanged.  This  term  is  important 
when  near  the  wall;  far  from  the  wall  its  contribution  is  almost  zero.  Its  role  is  to  transport 
turbulent  kinetic  energy  to  the  wall,  where  it  is  dissipated. 

7  Results 

The  initial  flow  field  was  laminar  with  small  random  fluctuations  in  two  dimensions  (zero 
divergence  determines  the  third  component).  The  calculation  was  allowed  to  run  until  the 
average  streamwise  velocity  profile  settled  down  to  steady,  fully  developed  turbulent  flow. 
At  the  grid  dimensions  mentioned  above  this  took  about  300  hours  of  Cray  2  time. 

The  computed  velocity  results  are  shown  in  Figure  4.  They  are  compared  with  the  hot 
wire  anemometer  data  of  Eckelman[21]  and  the  laser-Doppler  velocimetry  measurements  of 
Niederschultze[22].  If  the  data  were  plotted  with  semi-logarithmic  scales,  a  "logarithmic” 
layer  would  be  present  in  the  data  of  Eckelman,  but  not  in  the  other  two  sets  of  data.  This 
is  due  to  the  low  Reynolds  number  of  the  latter  sets  of  data  (roughly  a  value  of  2200), 
while  the  data  of  Eckelman  are  at  a  Reynolds  number  of  5600.  The  “logarithmic”  layer 
(see  Figure  1)  appears  at  a  Reynolds  number  of  about  4000. 

Turbulence  intensity  profiles  are  shown  in  Figure  5  and  Figure  6.  Another  second-order 
statistic  profile,  the  Reynolds  shear  stress,  is  shown  in  Figure  7.  The  streamwise  turbulence 
intensity  is  compared  to  Eckelman  and  to  Niederschultze;  note  the  maximum  in  the  center 
of  the  shear  layer  near  the  wall  at  a  distance  of  12  viscous  units.  The  normal  and  spanwise 
profiles  are  similar,  although  the  maxima  are  not  as  pronounced.  Good  agreement  exists 
between  the  computed  data  and  the  experimental  data,  although  the  maxima  are  slightly 
less  for  the  computed  data.  Hot  wire  anemometer  data  are  typically  high  when  measuring 
values  in  the  directions  perpendicular  to  the  bulk  velocity,  as  shown  in  the  figures.  The 
discrepancy  in  the  spanwise  maximum  can  not  be  explained. 

Skewness  and  flatness  profiles  are  shown  in  the  next  four  figures;  the  results  are  as 
expected.  Note  that  they  are  significantly  different  from  the  Gaussian  value  of  zero  (for 
skewness)  and  three  (for  flatness),  although  this  is  the  same  behavior  as  that  found  by 
Niederschultze.  The  values  of  the  flatness  near  the  wall  are  somewhat  large,  when  compared 
to  the  experimental  data.  Even  with  laser  velocimetry,  third-  and  fourth-order  statistics 
are  difficult  to  measure  close  to  the  wall.  However,  the  high  value  of  the  flatness  in  the 
streamwise  direction  suggests  that  the  simulation  does  not  match  nature  in  this  respect. 
Note  that  the  spanwise  skewness  and  the  normal  flatness  approach  their  Gaussian  values 
in  the  outer  flow. 

Figure  12  is  the  total  turbulent  kinetic  energy  budget,  as  described  by  Equation  31. 
Since  the  normal  and  spanwise  contributions  to  the  total  turbulent  kinetic  energy  budget 
are  very  small  when  compared  to  the  streamwise  component,  the  total  budget  is  very 
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similar  to  the  streamwise  budget.  The  terms  identified  on  the  figure  are  those  described 
in  Section  6.  Unfortunately,  there  are  no  experimental  data  to  compare  with;  however, 
the  values  are  very  similar  to  predictions  made  from  pure  shear  flows  and  extrapolated  to 
channels[2].  As  explained  in  the  introduction,  the  production  and  dissipation  of  turbulent 
kinetic  energy  (curves  V  and  T>)  are  such  that  the  transport  of  kinetic  energy  is  away 
from  the  viscous  wall  layer.  This  is  shown  in  the  Figure  12;  the  transport  terms  appear  as 
negative  quantities. 

Figure  13  shows  a  spanwise-normal  “slice”  of  the  instantaneous  velocities  in  the  span- 
wise  and  normal  directions.  Each  velocity  \rector  originates  from  one  of  the  grid  points  in 
the  solution  scheme  (the  grid  points  are  in  a  rectangular  array).  Along  the  walls  (i2=0 
or  300)  there  are  several  eddies,  which  appear  as  circular  masses  of  vectors,  as  well  as  the 
associated  in-  and  outflows  that  are  shown  schematically  in  Figure  2.  Note  the  outflow  at 
x2=0  near  x3=750  and  the  inflow  near  x3=300;  these  eddies  are  about  the  size  (50  units) 
that  would  produce  the  spacing  of  the  streaky  structures  discussed  in  Section  3.  In-  and 
outflows  do  not  appear  to  be  spaced  every  100  units  in  this  slice.  Future  work  will  involve 
averaging  techniques  to  examine  the  near-wall  turbulent  structure. 

Although  the  temperature  gradient  was  included  in  the  calculations,  there  are  no  sig¬ 
nificant  results  to  be  presented.  When  vector  plots  are  made  on  color-graphics  terminals, 
the  transport  of  heat  can  be  seen.  This  ability  to  track  other  transport  properties  of  the 
flow  field  may  be  of  interest  in  future  work. 

8  Conclusions 

A  direct  numerical  simulation  of  turbulent  flow  has  been  compared  with  experimental 
results.  The  computed  averaged  quantities  agreed  closely  with  those  values  obtained  in 
the  laboratory.  There  axe  discrepancies  in  the  maximum  values  of  the  second-  through 
fourth-order  statistics.  The  maximum  value  of  the  streamwise  turbulence  intensity  is 
about  ten  percent  lower  than  what  is  found  experimentally.  Since  the  experimental  data 
are  considered  very  accurate,  the  numerical  simulation  is  flawed  in  this  respect.  It  is 
possible  that  the  streamwise  length  of  the  computational  domain  is  not  large  enough  and 
structures  that  should  have  died  out  are  re-introduced  to  the  flow  field.  This  would  have 
the  effect  of  extending  the  correlation  length,  which  in  turn  would  lower  the  absolute 
values  of  the  zero-separation  correlation,  or  the  turbulence  intensity.  Since  the  computed 
turbulence  intensity  matches  the  experimental  values  in  the  normal  direction  (for  the  laser 
velocimetry  data),  the  streamwise  values  axe  even  more  suspect.  Currently,  efforts  are 
directed  to  increasing  the  computational  length  in  the  streamwise  direction  by  a  factor  of 
two,  while  maintaining  the  same  resolution.  These  results  should  be  available  for  the  next 
report. 
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Similar  arguments  hold  for  the  skewness  and  flatness,  although  the  experimental  data 
for  these  statistics  are  not  as  reliable  as  the  second-order  moments.  The  unusually  high 
flatness  at  the  wall  suggests  that  more  resolution  is  required  near  the  wall  to  sort  out  the 
contributions  of  various  sized  structures,  although  calculations  with  less  resolution  seem 
to  capture  the  smallest  relevant  scales. 

Future  direction  for  this  computer  simulation  involves  investigation  of  the  eddies  in 
the  viscous  wall  region.  This  work  will  center  around  some  type  of  averaging  technique. 
There  are  many  techniques  that  are  used  to  find  structures  in  wall  bounded  flow  fields[23]. 
Unfortunately,  many  of  these  techniques  are  akin  to  counting  wheels  on  vehicles  passing 
by  on  a  street  corner:  The  average  number  is  three,  but  not  one  tricycle  rode  past.  An 
averaging  technique  that  allows  counting  events  that  are  found  instantaneously  is  needed, 
instead  of  a  technique  that  results  in  an  average  structure  that  may  have  no  resemblance 
to  the  instantaneous  flow  field  at  all. 
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Expanded  View  of  the  Turbulence  Intensity 
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Dimensionless  Normal  Distance 
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Figure  13:  Instantaneous  Spanwise-Normal  Slice  of  the  Flow  Field 


Intentionally  left  blank. 
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Nomenclature 


The  general  system  of  nomenclature  followed  in  this  report  is  as  follows: 

1.  A  bold  face  symbol  represents  a  three-dimensional  vector  quantity. 

2.  A  subscript  of  1,  2,  or  3  represents  the  value  of  a  quantity  in  a  par¬ 
ticular  direction  where: 

(a)  Number  1  represents  the  streamwise  direction 

(b)  Number  2  represents  the  normal  direction 

(c)  Number  3  represents  the  spanwise  direction. 

3.  A  subscript  of  i  or  j  is  used  for  indicial  notation  of  vector  quantities 
and  second  degree  tensor  quantities  when  i  and  j  appear  together. 

4.  An  overbar  denotes  a  time  average. 

5.  A  tilde  denotes  a  fluctuating  quantity  (see  Equation  3). 

6.  A  plus  sign  superscript  denotes  a  dimensionless  quantity  (see  Equa¬ 
tion  8). 

7.  Upper-case  represents  the  mean  quantity. 

Lower-case  Symbols 

Symbol  Meaning 
h  Channel  half-height 

p  Hydrodynamic  pressure 

t  Time 

v  Velocity 

v*  Skin-friction  velocity 

x  Direction  in  the  computational  domain 


Upper-case  Symbols 


Symbol  Meaning _ 

A  Constant  for  temperature  equation 

Ao  The  van  Driest  parameter 

L  Number  of  Fourier  modes  in  the  streamwise  direction 

M  Number  of  Fourier  modes  in  the  spanwise  direction 

N  Degree  of  the  Chebyshev  polynomial  (normal  direction) 

Rij  Strain  rate 

Sij  Strain  rate 

T  Temperature 

T b+  Time  between  “bursts” 

U(x2)  Streamwise  component  of  the  mean  velocity 
Re  Reynolds  number 

Pr  Prandlt  number 

Greek  Symbols 


Symbol  Meaning 

k  The  von  Karman  parameter 

A  Length  of  the  computational  domain 

A+  Spanwise  spacing  of  the  streaky  structures 

u  Kinematic  viscosity 

ut  Eddy  viscosity 

p  Density 

r  Local  shear  stress 

TtJ  Reynolds  stress  tensor 

rw  Wall  shear  stress 

0  Transformed  temperature 

u >  Vorticity 

II  Modified  total  pressure 

II7'  Total  pressure 
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